{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import joblib\n",
    "import os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_01 = joblib.load('../../3.HA_model/1.PCA_T2/train_HA_feature/train_01.lz4')\n",
    "train_02 = joblib.load('../../3.HA_model/1.PCA_T2/train_HA_feature/train_02.lz4')\n",
    "train_03 = joblib.load('../../3.HA_model/1.PCA_T2/train_HA_feature/train_03.lz4')\n",
    "\n",
    "test_01 = joblib.load('../../3.HA_model/1.PCA_T2/test_HA_feature/test_01.lz4')\n",
    "test_02 = joblib.load('../../3.HA_model/1.PCA_T2/test_HA_feature/test_02.lz4')\n",
    "test_03 = joblib.load('../../3.HA_model/1.PCA_T2/test_HA_feature/test_03.lz4')\n",
    "test_04 = joblib.load('../../3.HA_model/1.PCA_T2/test_HA_feature/test_04.lz4')\n",
    "test_05 = joblib.load('../../3.HA_model/1.PCA_T2/test_HA_feature/test_05.lz4')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 312,
   "metadata": {},
   "outputs": [],
   "source": [
    "train = pd.concat([train_01, train_02, train_03], axis=0)\n",
    "y_train = train[['RULR']]\n",
    "# use_cols = ['PCA_T2__mean', 'PCA_T2__mean_diff', 'PCA_T2__delta', 'PCA_T2__delta_diff']\n",
    "use_cols = ['PCA_T2__mean', 'PCA_T2__mean_diff', 'PCA_T2__delta']\n",
    "x_train = train[use_cols]\n",
    "\n",
    "test = pd.concat([test_01, test_02, test_03, test_04, test_05], axis=0)\n",
    "# use_cols = ['PCA_T2__mean', 'PCA_T2__mean_diff', 'PCA_T2__delta', 'PCA_T2__delta_diff']\n",
    "use_cols = ['PCA_T2__mean', 'PCA_T2__mean_diff', 'PCA_T2__delta']\n",
    "x_test = test[use_cols]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 313,
   "metadata": {},
   "outputs": [],
   "source": [
    "from lightgbm.sklearn import LGBMRegressor\n",
    "\n",
    "reg = LGBMRegressor(boosting_type='gbdt', num_leaves=100, max_depth=50, learning_rate=0.1,\n",
    "              n_estimators=500, max_bin=100, objective='regression', \n",
    "              min_split_gain=0.01, min_child_weight=0.001, min_child_samples=1, subsample=1.0, \n",
    "              subsample_freq=1, colsample_bytree=0.9, reg_alpha=0.5, reg_lambda=0.5, \n",
    "              random_state=2018, n_jobs=12, silent=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 315,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LGBMRegressor(boosting_type='gbdt', class_weight=None, colsample_bytree=0.9,\n",
       "       learning_rate=0.1, max_bin=100, max_depth=50, min_child_samples=1,\n",
       "       min_child_weight=0.001, min_split_gain=0.01, n_estimators=500,\n",
       "       n_jobs=12, num_leaves=100, objective='regression',\n",
       "       random_state=2018, reg_alpha=0.5, reg_lambda=0.5, silent=False,\n",
       "       subsample=1.0, subsample_for_bin=200000, subsample_freq=1)"
      ]
     },
     "execution_count": 315,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "reg.fit(x_train, y_train.values.reshape(-1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 316,
   "metadata": {},
   "outputs": [],
   "source": [
    "RULR = reg.predict(x_test).reshape(-1,)\n",
    "CLR = 1 - RULR\n",
    "CL = test['CL'].values\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 317,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7f84cad1f5f8>"
      ]
     },
     "execution_count": 317,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xt83GWd6PHPN/dbc2mT3tJL0jYtlHIpxHIVWVBaXaUVWETUF15WjrsK58ihQHXP6uruWq1H3T0qLruy4q7QRcBaBakVRISVloTSKy1NL5SmlyRt0jTNbS7f88fvN8kkmWmTzH3m+3698srMb34zv+eXefKdZ57n+3seUVWMMcakr6xEF8AYY0xsWaA3xpg0Z4HeGGPSnAV6Y4xJcxbojTEmzVmgN8aYNGeB3hhj0pwFemOMSXMW6I0xJs3lJLoAAJWVlVpTU5PoYpg01djY2KaqVYk4ttVtE0ujrdtJEehrampoaGhIdDFMmhKRtxN1bKvbJpZGW7et68ZkLBF5RERaRGRH0LY1IrJbRLaJyC9EpDzosVUi0iQie0RkaWJKbczYWaBPgHVbmrl69QvUPvgMV69+gXVbmhNdpEz1E2DZsG0bgUWqehHwFrAKQEQWArcDF7jP+aGIZMevqCYTxCo2WKCPs3Vbmln19HaaO3pQoLmjh1VPb7dgnwCq+hJwcti236qq1737KjDDvb0cWKuqfap6AGgClsStsCbtxTI2WKCPszUb9tDj8Q3Z1uPxsWbDngSVyJzFp4HfuLergXeCHjvsbjMmYqrKN559M2axISkGYzPJkY6eMW03iSEiXwa8wM/G8dy7gLsAZs2aFeWSmXTg8ytvHu3ktYMnee3gSTYfaKetqy/kvtGIDRbo46y8KJf2bs+I7dPLCxNQGhOKiHwS+CBwgw6uzNMMzAzabYa7bQRVfRh4GKC+vt5W9jH0enxsO3zKDeonef3tdk73OT2E1eWFXDNvEi/uaaWjJzaxwQJ9HL30ViunejxkCfiD/v0Lc7NZuXRB4gpmBojIMuB+4D2q2h300HrgMRH5DjAdqAM2J6CIJgV09npoPNjO5oMnee3ASbYdPkW/zw9A3eQSPnTJdJbUTORdtROpdgN5oI8+uPsmWrHBAn2c7Gg+xV/9ZyPzp0zgzitr+P7vmzjS0cP08kJWLl3AisXW3RtvIvI4cB1QKSKHga/gZNnkAxtFBOBVVf2cqu4UkSeAXThdOp9XVV/oVzaZpqWzdyCobz7Yzu5jnahCTpawqLqMT15dw7tqJlI/u4KK4ryQrxGIAWs27Il6bJBkWDO2vr5e0/mikkMnurn5oVfIz8nm6b++iimlBYkuUkYRkUZVrU/EsdO9bmciVeXgiW43qDt97G+fcL78FeZmc+nsct5VM5ElNRO5ZFY5RXmxa0+Ptm5biz7GTnT1cee/b8bjU9be9S4L8sakmOEDp68dbKf1tDNwWlGUS33NRD5++WzeVTuRC6aXkpudfMmMFuhjqLvfy6cfbeBIRw8/+8vLmTd5QqKLZIw5h3MNnF49dxLvqnVa7HOrSsjKkgSX+Nws0MeI1+fn7se2sP1wBw99/DLqayYmukjGmBA6ez00vt3OawecFvvWdwYHTudPCT1wmmos0MeAqvLlX+zg+d0t/P2KRSy9YGqii2SMcQUGThsOtrP5wEneHMfAaaqxQB8D3/vdXv6r4R3uvn4eH79idqKLY0xGWLeleUTGyvJLpp9z4PR/3lAXl4HTRIrorETki8BfAgpsBz4FTAPWApOARuATqtofYTmTWnAFKyvMpaPHw19cNoN73zc/0UUzJiMMz0Fv7ujh3ife4G9/uYPOXqd/PVUGTmNh3IFeRKqBe4CFqtrj5hjfDnwA+K6qrhWRHwGfAR6KSmmT0PAK1uFeEHXFnEm4edjGmBgLNYeUX8HjU/7hw4tSauA0FiL9OMsBCkUkBygCjgLXA0+6jz8KrIjwGEktXAX7zsa3ElQiYzJPuPlgej0+Pnb5bOqmTMjYIA8RBHpVbQa+DRzCCfCncLpqOoKmeQ07w5+I3CUiDSLS0NraOt5iJJxNUmZM4oWbD8bmkHKMO9CLSAXOHN21OHN/FDNyEYewVPVhVa1X1fqqqoQs5xkVVsGMSbyVSxdQmDt0HRibQ2pQJF037wUOqGqrqnqAp4GrgXK3KwfOMsNfurAKZkzirVhczS2XDXYeZItwy2XVNoeUK5JAfwi4QkSKxBl1vAFnwqffA7e6+9wJ/DKyIia3FYur+cbNFzIh3/lsqy4v5Bs3X2gVzJg4WrelmacaB9uUPlWeamy2ldtckfTRb8IZdH0dJ7UyC2cO7geAe0WkCSfF8sdRKKcxxoRlK7edXUR59Kr6FZypXYPtJ4PW0gyVv7vq6e0A1qo3Jk6aLSnirDLjaoEYspaEMYnV6/GRF+bCJ0uKcFigj5ClV6YuEXlERFpEZEfQtokislFE9rq/K9ztIiL/LCJNIrJNRC5NXMlNgNfn5+7Ht+Dx+cnNHponb0kRgyzQR8jSK1PaTxiZEvwg8Lyq1gHPu/cB3o+zfGAdzsLfaXu1d6pQVb70i+1s3HWcr950AWtuvZjq8kIES4oYLj1n8ImjlUsXxGydRxNbqvqSiNQM27wcZ3lBcK7sfhEnwWA58FN3sfBXRaRcRKap6tH4lNYM983n9vBEw2HuuaGOO6+qAWxcLBxr0UcokL8b+NJo+bspb0pQ8D4GTHFvVwPvBO2X9ld9J7N/fWk/P/rDPj5+xSy++N66RBcn6Vmgj1Agfzew8q7l76YPt/U+5kWV0+Wq72T1ZONh/uHZN/nzi6bxdzctsskDR8ECfYQs6ybtHBeRaQDu7xZ3ezMwM2i/tL/qOxn9btdxHnhqG9fMq+Q7t11MdgZPVDYW1kcfIcu6STvrca7oXs3QK7vXA18QkbXA5cAp65+PveC1HiaV5HGqx8Oi6aX86BOXkZ+Tfe4XMIAF+ohNLy8MebGGZd0kPxF5HGfgtVJEDuNc/LcaeEJEPgO8Ddzm7v4szloLTUA3ziI7JoaGX4zY1tWPAH9RP5OSfAtdY2F/rQhZ1k3qUtWPhnnohhD7KvD52JbIBAvVLarAQy/usyU6x8gCfYQC2TX/Z90OTvd5qXbXqrSsG2MiY92i0WODscaYpGQXI0aPBfoIBfoRT/c5i2oFJjWz9EpjIrNy6YIRWTXWLTo+FugjZOmVxkTfui3NfPM3u/H5deBiRJvWYPysjz5C1o9oTHQNz7ZRBlvyFuTHx1r0EbJ+RGOi61vP7bZvyVFmgT5CtmasMZHz+5VX95/gS7/YzpFTvSH3sW/J4xdR142IlAP/BizC+Yb1aWAP8F9ADXAQuE1V2yMqZRKz9EpjxkdV2d58ivVvHOHX245yrLOXwtxsCnOzR7Towb4lRyLSPvp/Ap5T1VtFJA8oAr6EM5/3ahF5EGc+7wciPI4xJk00tXSxfusRfrX1CAfazpCbLbxnfhWrPnAe71s4hd/uPG4XIUbZuAO9iJQB1wKfBFDVfqBfRMLN552W1m1p5sGnttHr9QO2ZqwxoTR39PCrrUdY/8YRdh3tRASunDOJ/3HtHJYtmkp5Ud7AvoH/m8AcN9PtW3LEImnR1wKtwL+LyMVAI/A/CT+fd8rz+vwcOtnN3pYumtyfX287gsc3dCbbwMCRVUyTydq6+nh2+1HWv3GEhred3tuLZ5bztx9cyAcvmsbk0oKwz12x2NZ0iKZIAn0OcClwt6puEpF/YnDZNcCZH0REQs7nLSJ34SzJxqxZsyIoRvT1e/28feIMe1u62Hu8i70tp2lq6WJ/6xn6ff6B/arLC0cE+QAbODKZ6HSvhw07j7N+6xFeaWrD51fmTynhvhvn86GLpzN7UnGii5iRIgn0h4HDqrrJvf8kTqA/Hlhibdh83kOo6sPAwwD19fVjXtwhGno9Pva3nhkI5HuPd9HU2sXBtjN4/U6RRGBmRRF1k0t4z4Iq6iZPoG5yCXMnl1CSn8PVq1+w2StNRuv1+Hhhdwvr3zjCC3ta6Pf6mVFRyP+4dg43XTKd86aWJrqIGW/cgV5Vj4nIOyKyQFX34Mz4t8v9CTWfd8J093vZ1+IE9EArvanlNIdOduPGc7KzhNmTnIC+7IKp1E0pYd7kEuZUllCYF37ea5u90mQij8/PK01trN96hN/uPE5Xn5fKknzuWDKLmy6ZzuKZ5bbyUxKJNOvmbuBnbsbNfpw5urMIPZ/3uAUvPnC2gZnOXo/Td+52twSCenCLOzdbqK0s5oLpZSy/pJq6KSXUTZ5ATWXRuBYysIGj9CQiXwT+EidteDtO3Z4GrAUm4YxJfcJNQsgIfr/S8HY767c28+z2Y5w808+Eghw+cOFUll9SzRVzJtmKT0lKnGm2E6u+vl4bGhpCPjb8cmiAgpwsPnvtHKaWFQwMiu493sWxzsELLfJzsphbVeIG8hLmTZ5A3ZQSZk0sIjfbrhPLJCLSqKr1Y9i/GngZWKiqPSLyBIMLjzytqmtF5EfAVlV96Gyvdba6nYyGN6ruu3E+dVMmsH7rEX699QhHTvVSkJvFe8+fwk0XT+c9C6pspacEGm3dTvq5bkJNGtbr9fP/Xmgasq2yJI8bzpvMJTPLuWRWOQumTqC0IJf8nCz7CmnGIwcoFBEPzvUhR4HrgTvcxx8FvgqcNdCnkuGNquaOHu59YisK5GQ5ue4PvP883nv+FIpthaeUkvTv1mizV9q6+nl+dwvP7x469pudJRTnZVOcnzPwU5KfTVFeDiX5ORTnZ1OcF+6xoY+X5OdQkGsfHOlOVZtF5NvAIaAH+C1OV02Hqnrd3Q4DKd8/d7rXw77WMzS1dPHV9TtDruhUXpjL7++7jorivNAvYpJe0gf6sGuylhXwzD3v5ky/lzN9Prr6vJzp89Ld76Wrz8eZPq/72ODjwY+d6Ooe8tx+rz/E0UfKEoI+GNwPkGH3Swa2ZYf/AHE/fIrysu2DI8mISAWwHOdakQ7g58CyMTw/6VKHT3T1Dbn+Y1/ryO7OcE71eCzIp7ikD/ThslruX3YeFcV5UauAHp+f7j4fXf1euvu87gfH4AfEmb7QHyCB280dPUEfNF56PaP74BD3g6MoL3vgQyD4dvCHQqhvGUM+XPJzKMrNJssGxCL1XuCAqrYCiMjTwNVAuYjkuK36GUDI1WUSlTqsqhw51euOWZ1mX+tgYG/v9gzsV5SXzbzJJVw1dxJzJwfGsEr4+L9tCjmhmKUKp76kD/TxymrJzc6irCiLsqLcqLye1+en2+N+MLgfEgMfIO43iSGP9XsHvpWc6fNxrLN3yGPd/SMneQqnKOiDYfjtcB8gIx4L+raSgZkUh4ArRKQIp+vmBqAB+D1wK07mTcJSh70+P2+f7B5snbd0sddtpQfXk4qiXOomT2DZomnMc4P5vMklTC8rCPkt8v5l51mqcJpK+kAPqXk5dE52FqXZWZQWROeDw+dX95vF4LeIrj4v3e794A+JM8M+TLr6vLSc7nW+sQT2G8MHR2Fu9rBuqqAuqbwcivLH8AGSl01OFLKeRptyOx7uld5PAq8DXmALTgv9GWCtiPy9u+3HkR7rbOcRfEHfvhbnYr6mli4OtJ0ZckX2tLIC5k0u4bb6mc71H1VOQJ9Ukj+msliqcPpK+vRKExt+v9LjGfwg6O4f/BAI3A/1WHCX1cAHTZ+Xrn4vo61K+TlZzoeA2wUVbuB74ANk2IfLpv0n+N7v9tIXNK5SmJsddpm5saZXRtNYU4dzsoT5UybQ1eflnfbugb9plsDsScXMrSoZ0jqfW1XMhCg1JkzqSZv0ShMbWVkyEDgnR+H1VAMfHL6QHxaBbxEDHxbDBtE7ejwD4xyBbf4xtEFScSK5UKnDXr/y1vHTLF00lZsvrR4I6DWTiinItXx1Mz4W6E1UiAhFeTkU5eVQNWFsXQahqCp9Xv/At4bA2EZXn5dP/ftrIZ+TahPJhSuvz6/84I5L41wak84s0JukJCIU5GY7rdiSoY9Vh0u5TbHskLCpwyl2Hib52VwAJuWkyzq96XIeJvlZi96knHTJDkmX8zDJzwK9SUmpmHIbSrqch0luSZFeKSKtOFMaB6sE2hJQnGiz80i82apalYgDh6nbqSCV3+/RSJfzG1XdTopAH4qINCQq9zma7DxMKkr39zvdz284G4w1xpg0Z4HeGGPSXDIH+ocTXYAosfMwqSjd3+90P78hkraP3hhjTHQkc4veGGNMFFigN8aYNBf3QC8iy0Rkj4g0iciDIR7/nIhsF5E3RORlEVnobq8RkR53+xsi8qN4lz1EWcd1Lu5jF4nIn0Rkp7tPQXxLP6Sc431PPhb0frwhIn4RuST+Z2DG4lzvt7vPbSKyy62fjwVt9wW93+vjV+rRG0V9/m7QObwlIh1Bj90pInvdnzvjW/IYUtW4/QDZwD5gDpAHbAUWDtunNOj2TcBz7u0aYEc8yxvDc8kBtgEXu/cnAdmpdh7D9rkQ2Jfo98V+ovJ+1+EsrFLh3p8c9FhXos8h0vMbtv/dwCPu7YnAfvd3hXu7ItHnFI2feLfolwBNqrpfVftxlmRbHryDqnYG3S3GWYg+GUVyLjcC21R1q7vfCVUd/ZJP0RWt9+Sj7nNNcjvn+w18FviBqrYDqGpLnMsYidGcX7CPAo+7t5cCG1X1pHvuGxnDovDJLN6Bvhp4J+j+YXfbECLyeRHZB3wLuCfooVoR2SIifxCRd8e2qOcUybnMB1RENojI6yJyf8xLG16k70nARxj8hzHJazTv93xgvoi8IiKvikhwsCsQkQZ3+4pYF3YcRlWfAURkNlALvDDW56aapByMVdUfqOpc4AHgb9zNR4FZqroYuBd4TERKE1XG0QpzLjnANcDH3N8fFpEbElTEUQlzHgCIyOVAt6ruSEjhTLTl4HTfXIfT4v1XESl3H5utztQBdwDfE5G5iSliVNwOPJnAb9NxE9c8ehG5Eviqqi51768CmDRp0j/W1NTErRwmszQ2NrZpgiY1q6ysVKvbJlZGW7fjPU3xa0CdiNQCzTifqHfU1NT8oy0ObmJFRBI2e2RNTQ1Wt02sjLZun7PrRkQeEZEWEdkRtG2iiGx0U5A2ikiFu11E5J/dtKZtIjJk4UtV9QJfADYAbwJPqOrOsZyYMdESpm6vEZHdbv39RVCXBSKyyq3be0RkaWJKbczYjaaP/ieMHHl+EHheVeuA5937AO/H6durA+4CHhr+Yqr6rKrOV9W5qvoP4y24yWzrtjRz9eoXqH3wGa5e/QLrtjSP52V+wsi6vRFYpKoXAW8BqwDcawduBy5wn/NDEcnGmCiKUr0e4ZyBXlVfAk4O27wceNS9/SiwImj7T9XxKlAuItOiUlJjXOu2NPPgU9to7uhBgeaOHlY9vX3M/xSh6raq/tb95gnwKjDDvb0cWKuqfap6AGjCSeUzJirWbWlm1dPbI67XoYw362aKqh51bx8Dpri30zY9ySSPNRv20Ov1D9nW4/GxZsOeaB/q08Bv3NtWt03MnOr28LVf7aLHMzQBKFr1OuLBWFVVERlz6o6I3IXTvcOsWbMiLYbJIEc6esa0fTxE5MuAF/jZOJ5rdducVb/Xz+uH2nl5bxt/bGpj++EO/GGiaDTq9XgD/XERmaaqR92umcCVc83AzKD9ZrjbRlDVh3HnhK6vr0/Wq19NEppeXkhziMo/vbwwKq8vIp8EPgjcoIP5x1a3zbipKntbuvjj3jZe3tvKpgMn6e73kZ0lXDyjjC9cX8djm96mrat/xHOjUa/HG+jXA3cCq93fvwza/gURWQtcDpwK6uIxJipWLl3AA09toy+o+6YwN5uVSxdE/NruVaD3A+9R1e6gh9bjXKT3HWA6TsLB5ogPaNJWy+leXmlq449723ilqY3jnX0A1FYWc8ulM7imrpIr506itCAXgDmVxax6evuQ7pto1etzBnoReRznCrlKETkMfAUnwD8hIp/BWeH+Nnf3Z4EP4AxUdQOfiriExgyzYnE1z+8+zq+2Om2IbBFuuayaFYvH1mUepm6vAvKBjSIC8Kqqfk5Vd4rIE8AunC6dz2fCFZUmtHVbmlmzYQ9HOnqYXl7IyqULWHrBVDYdOMHLe9t4uamN3cdOA1BRlMtV8yp597xKrqmrZEZFUcjXDNTf4a871nodSlKsMFVfX692UYkZrXVbmkO26L9x84Uh/ylEpNG9bD/urG6nn0B2THDLO0tARPD5lbzsLOprKrimrpJ3z6vigumlZGVJTMoy2rod7ytjjYnYmg17hgR5GMxOiEbrx5iz+dZzu0dkx/gVSvKy+cHHLmVJzUQK85LrEgsL9CblxCPrxpjhTp7p57FNb3PkVG/Ix8/0eXnP/IRMqXROFuhNyol11o0xwZpaTvPjlw/y9OuH6fP6yc/JGvGNEpK7/lmgNyln5dIF3PfzNwj+X8vNlqhkJ5jME2pgdfkl03m5qY0fv3yAF/e0kpeTxc2Lq/n0NbXsOtIZs+yYWLFAb1LT8ByCxOcUmBQ0fGC1uaOHlU9u5ZvP7eboqV4qS/K5933z+djls5hUkg/A/CkTgNhkx8SKBXqTMrw+Px09Hv7x2TfxDgvsHr/aYKwZszUb9owYWPX4lNbTfXz7Ly7mQxdPIz9n5MDqisVjT+dNJAv0Ju5Ule5+H+3d/XR0e2jv7qe920NHdz/tZzzu9qBt7j6ne71nfV0bjDVjFa7OeP3KrZfNCPlYKrJAbyISaGUPBOQzIYL3sKB9qttDv2/kYFbAhIIcKoryqCjKpbwoj9rKYsqL8igvyqWiKI/v/e4t2rs9I56XzINhJjlNKysImUVTnWZ1yQK9AZxW9pl+37BA7dwO3B9rKzs3WygfFrAvLcob2FYRCN7Fg/uUFeaSm332SVXLCnNTbjDMJB+vz8+kkvwRgT4d65IF+gQKNdofjX4/j89PR7eHUz2xb2UHgnWgtR0cvIvzsnGnEYiqWF4qbtJb4H+uuaOHorxsuvt93Ly4mk0HTqZ1XbJAnyChRvtXPb0dGAxkZ2tlh+zf7u6n44yH033xb2XHW6oNhpnEG/4/193vIydLuHZ+Fd/5yCUJLl1sWaBPkFCj/T0eH/c/uY0fvtg0ELw9vvB5g8Gt7IqiPOYEtbIrip0gXV6YG5dWtjHJLtT/nDdDsrUs0CdIuNH+fp9/VK3s8sJccpKslW1MsjrQdibk1dSQGdlaFugTJNxl/NXlhfzLJxIy0aIxacXnV17Y3cJP/3SQP+5tC7tfJmRrRdQkFJEvishOEdkhIo+LSIGI1IrIJhFpEpH/EpG8aBU2naxcuoC8YS3ydBztT2Yi8oiItIjIjqBtE0Vko4jsdX9XuNtFRP7ZrdfbROTSxJXcnM2Jrj5+8Psmrv3W7/nsTxvYe7yLe983n6/ddAGFuUMvfsqU/7lxt+hFpBq4B1ioqj3uogy34yw88l1VXSsiPwI+AzwUldKmkRWLq3l2x1F+u/M4MP7FM0xEfgJ8H/hp0LYHgedVdbWIPOjefwB4P86qUnU4q6c95P42CTA8Y+2+G+czu7KY//jT2zyz7Sj9Pj9XzZ3E//ng+bz3/CkD3ZylhbkZma0VaddNDlAoIh6gCDgKXA/c4T7+KPBVLNCPsG5LMy/ubh2471PlqcZm6mdPzIiKlwxU9SURqRm2eTnOqlPg1N8XcQL9cuCn7hqyr4pIeWDd5PiU1gSEyli794mtKFCSn8NHl8zkE1fOZt7kCSOem6nZWuMO9KraLCLfBg4BPcBvgUagQ1UD+X2Hgcz7q47Cmg17RuSt2+IZSWFKUPA+Bkxxb1cD7wTtF6jbFujjLNTCH4pzId0rD15PSb4NPQ4XSddNBU4rpxboAH4OLBvD8+8C7gKYNWvWeIuRsmzxjOSnqioiY54XM9PrdrSpKvtau/jDW2384a3WsAt/dPZ4LMiHEclf5b3AAVVtBRCRp4GrgXIRyXFb9TOA5lBPVtWHgYfBWVczgnKkJFs8I2kdD3TJiMg0oMXd3gzMDNrP6nYMdfZ6+O8mJ7C/9FbbwP/KnKpiivOyOdM/cl12+98JL5JAfwi4QkSKcLpubgAagN8DtwJrgTuBX0ZayHS0cukC7n9y25Dum0zJAEhy63Hq7WqG1t/1wBdEZC3OIOwp65+PHr9f2d58ipfeauWlva28fqgDn1+ZkJ/DVfMm8dd/Npdr66qYObEo5OLc9r9zdpH00W8SkSeB1wEvsAWnFfMMsFZE/t7d9uNoFDTdWNZN4onI4zgDr5Uichj4Ck6Af0JEPgO8Ddzm7v4sTkZZE9ANfCruBU5R4eZ0ajndyx/d7piXm9o4eaYfgAury/ir98zl2vlVLJ5VPmL6DZvraOzESSJIrPr6em1oaEh0MeJq3ZbmkC36b9x8oVXYKBORRlVNyFVomVi3g4VqfedkCVNKCwa6YypL8ri2ropr51dxTV0lle5KTubcRlu3beQiQSzrxqQzv195+2Q3X/vVrpDzy7Se7uP+ZQu4tq6KhdNKycqy+ZdiyQJ9gljWjUkXvR4fe46dZtfRTnYd6WTX0U52H+0MOWAa4PH5+evr5sWxlJnNAn2CWNaNSUVtXX0DwTzwe39rF363B3hCfg7nTyvlL+pnsnBaKWs27KG1q2/E61g9jy8L9AliWTcmmfn8ysETZ4YE9TePdtJyejBoV5cXcv60Uj5w4TQWTivlgumlzKgoHDINdl5OlmXIJAEL9AliWTcmXs61kll3v5fdx04PBPU3j3ay++jpgeCckyXUTZnAu+uqOH/aBBZOL2XhtFLKi849X6FlyCQHC/QJYnPdmHgINS/M/U9uY+Ou44jArqOdHGg7QyD5rrQgh4XTS7l9idP1snB6KfMml5Cfk32Wo5xdps4vk0ws0MdYv9fP8c5ejnX2cvRUL8dO9XD0VC+Pbz5kWTcm5kKtqtTv8/PM9qPMqChk4bRSbrp4+kBQry4vtBXI0pAF+gj0enwcO+UG8M4eN5D3DvndFmIgqjgvm15P6IW4LevGRNPZ6tPLD1wfx5KYRLJAH8aZPm9QwO5xfncGB/Ie2rs9I55XVpjLtLICppYVsKi6lKmlhQP3A78nFORy9eoXLOvGxNzZVjIzmSMQDECnAAAWQUlEQVTjAr2q0tnrHRrAAwG9c7Br5XSvd8RzJxXnMbWsgOryAi6bXc60skKmlg4G8KllBRTlje5PunLpAstGMDF3343zuffnWwm+AN7qWeZJiUB/rqyBAFWlvdszMoAP61rpHnYhhwhUleQzrayA2spirppbOdgCLy1gWlkhk0vzKcgd/4DUcJaNYGIp8D8TaM0X5WXT0++zepahkj7Qh8sa+OPeVqomFAy0wAODnf3eoX3f2VnClAn5TC0r4PyppfzZgsnDulIKmTwhf8TESfFg2QgmFkLNL+P3K9/9yCVW3zJU0gf6cFkDT73eTG62OAG7tJCLZ5Sz7IKhAXxaWQGVJflk2zwaJoOEWoGp1+u3jK4MlvSB/mxZA3u+/n6bDMkY18kz/Ty++VDYFZgsoytzJX2gP1vWgAV5Y2D74VP85L8P8qttR+j3+snPyaLPOzJ91zK6MldEHdMiUi4iT4rIbhF5U0SuFJGJIrJRRPa6vysiOcbKpQsoHDYIalkDJtZE5IsislNEdojI4yJSICK1IrJJRJpE5L9E5NxzAMRIv9fPL99o5sM/fIUPff9lfrPjKB+pn8nGL17LN2+5yP5nzBCRtuj/CXhOVW91K30R8CXgeVVdLSIPAg8CD4z3ACsWV9Pw9kn+89VDgM0JY2JPRKqBe4CFqtojIk8At+OsMPVdVV0rIj8CPgM8FOvyBGedTSkt4KIZZbx+qIO2rj5qK4v5yocWcstlMygtyAWgbsoEwDK6zKBxB3oRKQOuBT4JoKr9QL+ILMdZng3gUeBFIgj067Y081Tj4BrMNieMiZMcoFBEPDgNmKPA9cAd7uOPAl8lxoF+eAbNsc5eju3qZeG0Uv7vbRfz7nmVIbswLaPLBIuk66YWaAX+XUS2iMi/iUgxMCVo0eRjwJRIChgq6yYwJ4wxsaCqzcC3gUM4Af4U0Ah0qGrgSrrDQMwj6erfjMygATjV4+E986tsnMqMSiSBPge4FHhIVRcDZ3C6aQaosyBtyEVpReQuEWkQkYbW1tZQuwC2EpOJP3dcaTlOY2Y6UAwsG8PzR1W3z0ZV+dXWIxzrtAwaE7lIAv1h4LCqbnLvP4kT+I+LyDQA93dLqCer6sOqWq+q9VVVVWEPEi5TwDIITAy9Fzigqq2q6gGeBq4GykUk0N05A2gO9eTR1u1wWjp7+dx/NnL341vIzQ7dYrf6b8Zi3H30qnpMRN4RkQWquge4Adjl/twJrHZ//zKSAq5cuoCVT27F4xv8YpCbLZZBYGLpEHCFiBQBPTh1uwH4PXArsJYo1G0YOtA6rayA6xZM5pntR+n1+PjSB85jUnE+f7Nuh82JZCISadbN3cDP3Iyb/cCncL4lPCEinwHeBm6L8BgjO39CdgYZEx2quklEngReB7zAFuBh4BlgrYj8vbvtx5EcZ/hA65FTvTy2+RC1lcX8+M565lSVAM40HpZBYyIRUaBX1TeA+hAP3RDJ6wZbs2EPHv/QyO7xq13ObWJKVb8CfGXY5v3AkmgdI1SiAUCfxzcQ5MEyaEzk4j+T1xjZYKxJV+Hq8NEwUxgYM15JH+htMNakK6vbJl6SPtDbFAgmXVndNvGS9JOa2QIdJl1Z3TbxkvSBHmwwyqQvq9smHpK+68YYY0xkRDXxSeki0oqTcx8PlUBbnI4VK3YOYzNbVcd+iWoUjLJuJ+r9TGQ9ysRjx+K4o6rbSRHo40lEGlQ1VO5/yrBzSC+J+lsk8j3IxGMn8pyt68YYY9KcBXpjjElzmRjoH050AaLAziG9JOpvkcj3IBOPnbBzzrg+emOMyTSZ2KI3xpiMktKBXkSWicgeEWlyFyIf/vjnRGS7iLwhIi+LyEJ3+xJ32xsislVEPhz0nINBz2lI1nMIenyWiHSJyH2jfc0UOo+4vhfRNtr3QURuEREVkfqgbavc5+0RkaXxOraI1IhIT9D/x4+ifWwR+aSItAYd4y+DHrtTRPa6P3fG8bi+oO3ro33O7j63icguEdkpIo8FbR/3OY+aqqbkD5AN7APmAHnAVmDhsH1Kg27fBDzn3i4CctzbgVWwAvcPApXJfg5B254Efg7cN9rXTIXziPd7kYi/ibvfBOAl4FWg3t220N0/H2c5w31AdpyOXQPsiHFd+CTw/RDPnYgzFfREoMK9XRHr47qPdcX4nOtw1jCocO9PjvScx/KTyi36JUCTqu5X1X6cVX+WB++gqp1Bd4txlyxR1W4dXOS5gMQtZTLucwAQkRXAAWDnWF4zBmJxHqlutO/D14FvAsFzEy8H1qpqn6oeAJoY2zz4kRw7UpHUv6XARlU9qartwEZGv1ZvIur9WI79WeAH7nmhqoElViM551FL5UBfDbwTdP+wu20IEfm8iOwDvgXcE7T9chHZCWwHPhcU+BX4rYg0ishdMSu9Y9znICIlwAPA343nNaMsFucB8X0vou2cfxMRuRSYqarPjPW5MTw2QK2IbBGRP4jIu8dw3FEd23WLiGwTkSdFZOYYnxvt4wIUiLOg+6tuw2MsRnPs+cB8EXnFPcayMTw3Yqkc6EdFVX+gqnNxgsnfBG3fpKoXAO8CVolIgfvQNap6KfB+4PMicm3cCz1MmHP4KvBdVe1KWMHGaBznkXTvRbSISBbwHeB/J9mxjwKzVHUxcC/wmIiURrkIvwJqVPUinBbso1F+/fEcd7Y6V63eAXxPROZG+dg5ON031wEfBf5VRMqjfIywUjnQNwPBn8gz3G3hrAVGfFKr6ptAF7DIvd/s/m4BfkEUl44LIZJzuBz4logcBP4X8CUR+cI4XjMaYnEe8X4vou1cf5MJOHXuRffcrwDWu4Oikb6H4z622110AkBVG3H6nudH8dio6glV7XPv/htw2WifG6PjBte1/cCLwOJRHne05T4MrFdVj9sd9xZO4I/L/2tS5NFXVlZqTU1Nooth0lRjY2ObxnlSMxHJwflnvgHnH/c14A5VDTkOISIv4gxEN4jIBcBjOB9s04HngTpVHbnAbPSPXQWcVFWfiMwB/ghcqKono3VsEZmmqkfd2x8GHlDVK0RkItAIXOru+jpw2WiOHeFxK4BuVe0TkUrgT8ByVd0VxXNeBnxUVe90j7EFuASne3Jc5zwWSTEffU1NDQ0NKZc9Z1KEiMRrZtQBqup1v5lswMnKeERVd4rI14AGVQ2bwufu9wSwC/ACnx9tkI/02MC1wNdExAP4ccavRh10Rnnse0TkJvfcTuJkw6CqJ0Xk6ziBEuBroz12JMcFzgf+RUT8OL0cq0cb5Mdw7A3AjSKyC/ABKwPfnMZ7zmORFC36+vp6tUBvxmLdluZRr8wkIo1qM2WaDJYULXpjxmLdlmZWPb2dHo/TyG3u6GHV09sBbLUmY0JI5cFYk6G+9dzugSAf0OPxsWbDngSVyJjkZi16k3R6PT4Ot/dwuL2bw+09NHf0DLnferov5POOdPTEuaTGpAYL9Cbuuvu9NLe7wbtjMIAfbu+hub2btq7+IfvnZgvTywuZUVHI9Qsm8+yOo5zu9Y543enlhfE6BWNSigV6E3Vn+rxuS3wwgA+0ztt7OHFmaCDPy86iusIJ5AsXTqG6vJAZFUXMqHB+V03IJztLBva/cu6kIX30AIW52axcuiBu52hMKrFAn6HGkrUy3Olej9OdcnJYa9xtnbd3e4bsn5eTNRC0F1WXuYHcuT+zopDKknyyggL5uQTKOd7yG5NpLNBnoHNlrXT2ekIE8cHbp3qGBvKC3KyBFvhFM8qCWuOFVFcUUlk8tkA+GisWV1tgN2aULI8+A129+gWaQwxc5mYJBXnZI/q/i/KynaA9rEslEMgnFechEt1AHk2WR28ynbXoM1CoIA/g8St3LK4eEsRnVBRRUZSb1IHcGHN2FugzSJ/Xx3/86W2yBPwhvshVlxfyd8sXxb9gxpiYskCfAXx+5Rdbmvnuxrdo7uhhwZQJHDxxhj6vf2Afy1oxJn2dM9CLyCPAB4EWVV0UtP1u4PM4E/Q8o6r3u9tXAZ9xt9+jqhtiUXAz0vBMmvtunE9JQS5rNuzmreNdXDSjjG/dehFXz6uMKOvGGJNazjkY6y720AX8NBDoReTPgC8Df+5O7TlZVVvEWfD5cQanV/0dMP9cM+/ZYGzkhmfSAANdNHMqi7lv6QLev2hqRva122CsyXTnbNGr6ksiUjNs81/hTOXZ5+4TWP9wYK1L4ICIBNa6/FPUSmxCWrNhz4j5X/wK5YW5bPjiteRm27RGxmSq8f73zwfeLSKb3HUl3+VuH/X6hyJyl7tGY0Nra+s4i2ECws3zcqrHY0HemAw33giQA0zEWYJsJfCEjLFPQFUfVtV6Va2vqorr4j9pKdw8Lzb/izFmvIH+MPC0OjbjrERTSWLWKzXAyqULKMzNHrLNMmmMMTD+QL8O+DMAEZkP5AFtwHrgdhHJF5FanMVvN0ejoObsViyu5pbLqgfe0CzglstsmgBjzCgCvYg8jjOYukBEDovIZ4BHgDkisgNYC9zptu53AoG1Lp9jjGtdmvFbt6WZpxqbCWTG+4GnGptZt8W+UBmT6WyumzQRbv6a6vJCXnnw+gSUKHlYeqXJdJaOkSbCZd3YqkvGGAv0acKybowx4VigTxOWdWOMCccCfZqwrBtjTDgW6NOEZd0YY8KxQJ8GVJV/fPbNEXPd9Hh8rNmwJ0GlMsYkC5uPPgX5/MqbRzvZdOAkmw+c4LWD7Zw80x9yX8u6McZYoE8B/V4/25tPsfnASTYdOEHjwXZO9znrus6aWMT1503md7uO0zFs0W6wrBtjjAX6pNTT72PLO+1sPnCSzQdO8vqhdno9Tu/7vMklfOiS6VxeO5EltROZVuYE8lDz0VvWjTEGLNDHVbhVnU73emh8ezCwbz3cgceniMDCaaV8dMksLq+dSH3NRCpL8kO+diC7xlaNMsYMZ1MgxEmoFndOljCtrIDmjh786ty/cEYZS2oncnntRC6bPZGywtwEljo92BQIJtNZiz4OWjp7+btf7RyRFeP1K8dP9/GF6+u4vHYii2eVU5Rnb4kxJrosqkSZqnK4vWegG2bzwZMcaDsTdn+P18+975sfxxIaYzKNBfoIqSr72844GTH7T7D5wEmOnOoFoKwwl3fVTOSOJbN4+KX9tHb1jXi+ZcUYY2LtnIFeRB4BPgi0qOqiYY/9b+DbQJWqtrnLCf4T8AGgG/ikqr4e/WInjt+v7D52ms0HTrD5oNNqb+tyctgrS/K5vHYin5vjZMTMnzyBrCxnhcWqCfmWFWOMSYjRtOh/Anwf+GnwRhGZCdwIHAra/H6cVaXqgMuBh9zfKSFUVsyfXzSNnUc6ncDudsd09jo57NXlhVxbV8USN9WxtrKYcEvnWlaMMSZRRpV1IyI1wK+DW/Qi8iTwdeCXQL3bov8X4EVVfdzdZw9wnaoePdvrJ0PWTaismCyBnOws+r1ODvucyuKBoL6kdiIzKooSVVwzBpZ1YzLduProRWQ50KyqW4e1YKuBd4LuH3a3nTXQJ5LX52db8yn+9pc7RmTF+BVys4Tv3LGYJbUTmTyhIEGlNMaY8RtzoBeRIuBLON024yYidwF3AcyaNSuSlxoTv1/ZdbSTP+07wX/va+O1g+10udMJhNLd7+ODF02PW/mMMSbaxtOinwvUAoHW/AzgdRFZAjQDM4P2neFuG0FVHwYeBqfrZhzlGBVVpamli//ed4I/7TvBqwdO0NHtzAkzp7KY5ZdM56q5lXz917s41tk74vmWFWOMSXVjDvSquh2YHLgvIgcZ7KNfD3xBRNbiDMKeOlf//FiFm0YgqHwcOtntttidnzY3rbG6vJD3nT+Fq+ZN4so5lUwtG+yK8fj8lhVjjElLo0mvfBy4DqgUkcPAV1T1x2F2fxYntbIJJ73yU1EqJzBywLS5o4dVT2+no7uf0sLcgVZ7szs1b9WEfK6eN4mr5k7iqrmVzJwYfvDUsmKMMekqpea6uXr1CwNBPJTyolyunOME9ivnVjK3Kny6o8kclnVjMl1KXRl7tkU0nrnnGs6fWjpwgZIxxhhHSi0lGG5gtLq8kAuml1mQN8aYEFIq0K9cuoDC3Owh22zA1Bhjzi6lAv2KxdXccln1QKGzgFsuq7YBU2OMOYuUCvTrtjTzVGMzfve+H3iqsZl1W0Km6htjjCHFAv2aDXtGTFPQ4/GxZsOeBJXIGGOSX0oF+nBZN2fLxjHGmEyXUoE+XNaNTVNgjDHhpVSgv+/G+QzPoLSsG2OMObuUuGAqML9N4KrYwtwsej1+m6bAGGNGIekDfagFQVThux+5xAK8McaMQtJ33YTKtOn1+i3TxhhjRinpA71l2hhjTGSSPtBbpo0xxkTmnIFeRB4RkRYR2RG0bY2I7BaRbSLyCxEpD3pslYg0icgeEVkaaQFtfhtjjInMaFr0PwGWDdu2EVikqhcBbwGrAERkIXA7cIH7nB+KSDYRWLG4mm/cfCHV5YUIzkyV37j5QhuINcaYUTpn1o2qviQiNcO2/Tbo7qvAre7t5cBaVe0DDohIE7AE+FMkhVyx2CYuM8aY8YpGH/2ngd+4t6uBd4IeO+xuM8YYkyAR5dGLyJcBL/CzcTz3LuAu926XiCQiX7ISaEvAcccqVcoJyVnW2YkugDGJNO5ALyKfBD4I3KCDC882AzODdpvhbhtBVR8GHh7v8aNBRBpSYS3RVCknpFZZjckU4+q6EZFlwP3ATaraHfTQeuB2EckXkVqgDtgceTGNMcaM1zlb9CLyOHAdUCkih4Gv4GTZ5AMbRQTgVVX9nKruFJEngF04XTqfV1Vf6Fc2xhgTDzLY65J5ROQutwspqaVKOSG1ympMpsjoQG+MMZkg6adAMMYYE5m0DPQissydgqFJRB4Ms89tIrJLRHaKyGNB230i8ob7sz7RZRWR7waV5y0R6Qh67E4R2ev+3JnE5Yzr39QYM1Tadd24Uy68BbwP54Kt14CPququoH3qgCeA61W1XUQmq2qL+1iXqpYkS1mH7X83sFhVPy0iE4EGoB5QoBG4TFXbk6mc7v24/U2NMSOlY4t+CdCkqvtVtR9YizM1Q7DPAj8IBMVAkE+A0ZQ12EeBx93bS4GNqnrSPY+NjJyTKBnKaYxJsHQM9KOZhmE+MF9EXhGRV93rAgIKRKTB3b4iCcoKgIjMBmqBF8b63CiIpJwQ37+pMWaYpF9KMEZycC7mug7n6t2XRORCVe0AZqtqs4jMAV4Qke2qui+BZQ24HXgyBa5LCFXOZP2bGpMR0rFFP5ppGA4D61XVo6oHcPqf6wBUtdn9vR94EVic4LIG3M7Q7pCxPDdSkZQz3n9TY8ww6RjoXwPqRKRWRPJwAs/wTI91OK15RKQSpytnv4hUiEh+0Parca7yTWRZEZHzgAqGTve8AbjRLXMFcKO7LanKmYC/qTFmmLTrulFVr4h8ASfoZQOPuFMzfA1oUNX1DAbJXYAPWKmqJ0TkKuBfRMSP8yG4OlxmSRzLCk5gXRs0eRyqelJEvo4ThAG+pqonk62cwPnE8W9qjBkp7dIrjTHGDJWOXTfGGGOCWKA3xpg0Z4HeGGPSnAV6Y4xJcxbojTEmzVmgN8aYNGeB3hhj0pwFemOMSXP/H+QrwJylM6vkAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 5 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "plt.subplot(3,2,1)\n",
    "plt.plot(CLR[:8],CL[:8])\n",
    "plt.scatter(CLR[:8],CL[:8])\n",
    "\n",
    "plt.subplot(3,2,2)\n",
    "plt.plot(CLR[8:15],CL[8:15])\n",
    "plt.scatter(CLR[8:15],CL[8:15])\n",
    "\n",
    "plt.subplot(3,2,3)\n",
    "plt.plot(CLR[15:25],CL[15:25])\n",
    "plt.scatter(CLR[15:25],CL[15:25])\n",
    "\n",
    "plt.subplot(3,2,4)\n",
    "plt.plot(CLR[25:35],CL[25:35])\n",
    "plt.scatter(CLR[25:35],CL[25:35])\n",
    "\n",
    "plt.subplot(3,2,5)\n",
    "plt.plot(CLR[35:45],CL[35:45])\n",
    "plt.scatter(CLR[35:45],CL[35:45])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 318,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import curve_fit\n",
    "def linear_fit(x, a):\n",
    "    return a*x\n",
    "def plot_fit(x, y):\n",
    "    x_plot = np.arange(0,1,0.01)\n",
    "    x = np.hstack((x,[0]))\n",
    "    y = np.hstack((y,[0]))\n",
    "    popt, pcov = curve_fit(linear_fit, x, y)\n",
    "    y_plot = linear_fit(x_plot,*popt) #拟合y值\n",
    "    plt.plot(x, y)\n",
    "    plt.scatter(x, y)\n",
    "    plt.plot(x_plot, y_plot)\n",
    "    return linear_fit(0.98, *popt)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 319,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl8VPW5+PHPk8nKmrBDIIR9FwJB3FkVcAO1Vay4YnFtb5dL1ba3Wm/7k6ut1mqrUkBBraKoiFb2VUU0CRHCFiCELeyEhIRsk5nv749zgkNIyCSZyUyS5/16zSsz55w588zhnIcz31WMMSillGr4QgIdgFJKqbqhCV8ppRoJTfhKKdVIaMJXSqlGQhO+Uko1EprwlVKqkdCEr5RSjYQmfKWUaiQ04SulVCMRGugAANq0aWPi4+MDHYZqoFJSUk4aY9oG4rP13Fb+VN1zOygSfnx8PMnJyYEOQzVQIrI/UJ+t57byp+qe21qko+q/rE2w6lnQcaFUQ+Isgs0LYOtHPttlUNzhK1VtzkLY+jEkzYbDmyCsKQy5C1r3CHRkStVO9l5IfhNS34HCbOg5Dgbe5pNda8JX9cupDEiea10MRTnQpjdM+D8YcidEtgx0dErVjKsUdi2F5DmQsRrEAX2vh8Rp0G2kzz5GE74Kfm4X7FoGSf+yLoaQUOh7Awx/EOKvBpFAR6hUzZw5ApvmQco8yDsMLWJh1G9h6D3QoqPPP04Tvgpe+cd/uBhyD0LzjjDqKRh6r18uBqXqhNsNmWutX6o7vwDjgh5j4Ya/QK/x4PBfWtaEr4KLMXDgG6tsfvticDutn7QTnoPeE/16MSjlVwXZ8P27VqLP3gtNWsPlj0Hi/dCqe52EoFePCg5FZ2DLAutiOL4dIlpaRTbDp0GbXoGOTqmaMQYOJUHSHNj2CbiKIe5yq9im/80QGlGn4WjCV4F1bJt1MWxZACX50HEw3PyK1SohvGmgo1ONzKLULF5Yls7hnEI6RUcxY3wfJifEVn9HxXmw5QOrtc2xNAhvDkPvhsQHoP0A3wfuJU34qu6VlsCOxVaiP7ABHBFWgh8+DWKHaSWsCohFqVk89XEahU4XAFk5hTz1cRqA90n/6FbrV2rZDUyHQXDj32DQjyGimb9C95omfFV3cg5CyltWRezZExATD9c+Cwl3Q5NWgY5ONXL/t3TnuWRfptDp4oVl6RdP+M4i2P6p1aTy4Lf2DcytVpPKzolBdQOjCV/5l9sNe1dbd/O7llrLeo237uZ7jIUQ7eytAutEXjGvrN7NkdyiCtcfzims+I2nMiDlTUh91+og1ao7XPcnqwNgkN7AaMJX/lG+RULTtnDVL2HYfRAdF+joVCO3KDWL577YwbG84iq37RQd9cOLqjpIBfkNjCZ85TvGWOPaJM22xv8oa5Ew+nfQ76Y6b5GgVEVeXb2bl1bsxlVu7KUJA9uzLv3kecU6UWEOZozvA2cOw6b5ddZByl804avaKymwEnzSbDjyPYQ3g4S7rGaVAWyRoFQZYwwb92bz2roM1u86UeE2aYfO8Nytg8610oltGcH/DT3Nlem/gcVL6rSDlL/Uv4hV8Di5x/pp+/27UJQLbfvBDX+FS+6AiOaBjq5KIjIXuBE4bowZaC9rBSwA4oF9wO3GmNMiIsDLwPVAAXCfMWZTIOJW3nO7Dcu3H+P1dRl8fzCHNs3CK902K6eQyQmxTO4TZY3VlPImfLMXolrVeQcpf9GEr6qnrAwz6V+wdy2EhFnFNcOnQdcrg6pFghfeAl4F5nssexJYZYyZKSJP2q+fACYCvezHCOA1+68KQiWlbhZ9n8Ub6zLIOHGWuFZN+NPkgfxoWGcG/GHZBcU5YEgM2QMfP/RDB6kul1lDefSf1GCKIzXhK+/kHbXLMN+CM1nQojOM/r1Vhtm8faCjqxFjzHoRiS+3eBIwyn4+D1iLlfAnAfONMQbYKCLRItLRGHOkbqJV3sgvLuX97w4w+8tMjp4pol/HFvz9zgSuH9iBUIdVoeqZ7JtSyCTHBqY6VtI/ZD/sbA4JU60bmAZYHKkJX1XOGNj3lVU2v/NzcJdaZZjXv1BvyzC90N4jiR8Fyv43iwUOemx3yF6mCT8InMov5q0N+5i3YR9nikq5rHsr/u9Hl3BNrzZIuV+dMU3CaFeYwVTHSiY7vqa5FLLd3ZXnQx/mN7/+fb0ojqypBnnFqloqyrVm2kmaDSfTITIaRjxsdQtvRBOMGGOMiFR7Gi0RmQ5MB4iL0yao/nQwu4B/fbmXD5IPUlzq5rr+7Xl4ZA8S4mIu3NhZRPKSN/lX6ZskRqRTbML43H0Z75SOY2tIL164YUiDTvagCV95Oppmj2vzATjPQqehMOmfMOAWCG8S6OjqyrGyohoR6Qgct5dnAV08tutsL7uAMWYWMAsgMTFR5130gx1HzvDGugw+23KEEIFbEmKZfk0PerarYPgCjxmkEguz2UsH/td5FwtdI8nF2j46PKxmY+bUM5rwG7vSYqtbeNJsq1t4aCQM+pHVkSR2aKCjC4TFwL3ATPvvpx7LHxeR97Eqa3O1/L5uGWNI2nea19buYU36CZqEO7j/inimXd2Nji2jzt+4gg5SWR3GMiM3kW/c/THlpvPOLXTW4TcJHE34jdXp/Vazs01vQ8FJaNUDrvszDPlJ0HYL9zUReQ+rgraNiBwCnsZK9B+IyDRgP3C7vfkXWE0y92A1y7y/zgNupNxuw6qdx3lt7R42HcihVdNwfn1tb+6+vCvRTco1s7zIDFK3/3MHWe6Kh0k4rzdtA6YJvzFxu2HPSutufvdyqwlln+utFgndRgV9t3BfM8bcWcmqsRVsa4DH/BuR8lRS6mbx5sO8sS6D3cfz6RwTxbOTBvDjYV2ICnf8sOFFZpD69OxAnl+RweGlm7hY2dqM8X38/XWCgib8xuDsKUh927ogcvZD03ZwzQwYdi+07Bzo6JQ6z9niUt5POsicL/dyOLeIvh2a8/KUIdwwqOO5ppVAxTNIXfE4DLsfWnWzhjtelHbBCJjlRUc1jvJ78CLha2/EesoYOJRs3c1v+xhcJdD1Khj3NPS9CUIr73GoVCBkny1h3oZ9zPtmHzkFTi7t1oo/3zKIUX3a/tC0stIZpKwOUovSTvLCrHQO52wnRKSCDlbniwpz8MzNDa+9fWW8ucN/C+2NWH+UnIW0D60L4ugWe6ade61im3b9Ah2dUhc4dLqA2V9msiDpIIVOF+P6teeRUd0Z1tWjLsmLGaTKT2BysWQvULsZreqpKhO+9kasJ07ssse1+TcUn4H2A+GGF+GS2xt822JVP6UfzeONdRks3nwYgElDYnl4ZHd6tfc4X8tPgdl+ENz4Egy6/YIZpF5Yll5l8Q1AbHQUXz85xqffpb6oaRm+9kYMBi4n7PyPVWyz70trXJv+k6xRKuMuq2/j2qhGInlfNq+tzWDVzuNEhTm453KraWVsWUuZczNIzYWDGz1mkHoAOg+v9LyudKISD+eGO26kal1pq70Rq+aziZHLnDlsNTtLeQvyj0LLOBj7B0i4B5q19VncSvmKMYY16cd5bW0GSftOE9MkjF+O6809l3clpqldn5S910ry52aQql5T4U7RUWRVkPQdIriNaZRFOOXVNOFrb0Qv+WRiZLAqqzLX2ePafAHGDT3HwfC/Qa/rIMRR9T6UqmNOl5vPtxzm9bV7ST+WR2x0FE/f1J87hnehSXio1UFqx+c+mUFqxvg+511rYN3RP3froEad5D3VNOFrb0QvVVSu6NXEyOc2zoHN71nlmKd2N6ixuVXDVVjiYkHSAf71ZSZZOYX0bt+MF28fzE2DOxHmCLF+pW7wmEGqeadazyBVdj359Nd0A+NNs0ztjVhDO46cqfAnJnhR3nj4e+uuJ20hOAsgNhEmv26NaxMW6Ydolaq902dLmP/NfuZ9s4/ssyUkdo3h2UkDGN2nHSEYq4NU0hxIL5tBaow1+mrvCT4ZfXVyQqwm+IvwppWO9kaspsM5hfx1+S4+Tj2EQIU9/KKbhF240FlktS1OnmO1NQ5rAoN+bFVWdRri77CVqrHDOYXM/jKT95MOUFDiYmzfdjw8qgfD41tZHaQ2vmo1qczO0F+pAaQ9bX0ot9DJP9fu4c2v9wEw/eruzN+QSWHphSm/2LOYJzvTrqx6x6qsat0TJsyEwXdCVHQdRa9U9e05nsfr6/ayKDULA0wa3ImHRvagT/tm1k3LJ3Nh68ceM0g9Cf1u1l+pAaIJ3weKS128/c1+Xl2zh9xCJ7cMieVX1/Wmc0wT3li/t8L3FDlLrZ+1SXOs8W0kxKqsGv6gVVmlTSpVENt04DSvrc1gxfZjRIaFMPWyrjx4dTc6N3FD2gdWoq+kg5QKHE34teB2Gz7bcpgXlqVz6HQhV/dqw5MT+zKgU8tK39OaXO5wrOUnoavgvZPQrAOMfMIa16ZFpzqMXqnqMcawdtcJXlubwXeZ2bSMCuPnY3tx7+VdaX12D3z9e6s3bEme3UHqb1aRZEQFY9SrgNCEX0Mb9pzk/y3ZwdasM/Tr2IL5Dwzimt6VtYE3JEo6U0NXMjHkOyKklA2u/nSe8lfoewM4KijPVypIlLrc/CftCK+tzWDn0Tw6tozkf27sz5SEtjTN+AIW/KxcB6lp0DlRf6UGIU341bTz6BlmLtnJ2vQTxEZH8dIdg5k0OJaQkApO7uJ8fuJYxd2OFfQLOcAZE8W/XWN5xzWODBPLvgE31P0XUMpLRU4XHyYf5I31ezl0upCe7Zrxlx8P5uYuxYR/PxdefadGHaRU4GjC99LhnEJeXLGLjzYdonlEKL+9vi/3XB5PZFgFHZ6O77DK5je/z/8Ly2ObuytPOH/KYtflFKKVVSq45RY4eXvjPt78eh+nzpaQEBfN09f3ZqwjlZCUn8FntesgpQJHE34VcgudvLY2gze/zsQY+OnV3Xl0VI8LZ9opLYGdn1uJfv9X4AjnYMfx/FfGMDaZXljj8/1Af+yqYHM0t4g5X+3l398e4GyJi1F92vLz4c1IOPEpsnz+BTNI1bSDlAocTfiVuFjLm/PkHrLGtNk0H/KPQXRXGPcMJNzDlFc2k2Uq7mBlsIZd0E4iKtD2HM9n1voMPknNwuU23HxJB37Z4whdM/8OH50/gxS9xvukg5QKDP2XK6eiljdPTOjLwNiWnht59Bj8whrnptd1VpPKnmPPjWtTVW/aGo2po5SPfH8wh9fW7mH59mOEO0J4YGhLHmn5LdHbfw/p9gxS2kGqQdGE72HDnpM8t2QnaVm5Fbe8Kci2xptPnuMxpdrPrTbGMV0v2F9lo/eVqdaYOkr5gDGGL3ef5LW1GXyz9xQtIh38ObGAW13LiNy++IIZpAiNCHTIyoc04XN+y5tOLSN58fbBTB7i0fIma5N1N791IZQWQZcRMPJJGDD5ohdERaP3lefNGN5K1Vapy82SrUd5fV0G2w6foVtzN+8M3s7lpxfjSNsK4c0gYao1M5p2kGqwGnXCP5JbyIvLd7GwopY3zkLY/LE1HPHhTRDW1BrqYPg06DDIq/17jt5X2Z1+p7JJH5TygyKni4Uph5i1fi8HsgsY1+oEL/feQI+j/0HSPWeQ+rHOjNYINMqEf6bIankz9yur5c2DV3XjsdE9rZY3pzJ+GNemKAfa9oWJL8DgOyCy8h60lSkbva/8uPigs+8o/8ktdPLOxv28+fU+8vLzeKjtVu6PXUXMqVTI8m4GKdXwNKqEX1zq4t2NB3hl9W5OFzi5JSGWX13bmy4tw2H3MutuPmM1hIRC3xutStj4q3xyQehY3aouHDtTxNyvMnn32wO0Ksni2TYbuNaxkrC801bF63V/giF3aQepRqpRJHy32/B52hFeWLaTg9mFXNXTGvNmYIsiSP0HJL8FZw5ZkzCM/p3Vxrh5B5/HoWN1K3/JPHmWWeszWJRygGtI5sOWX9FPkiBfO0ipHzT4hL8h4yQzl+xkyyG75c39A7kmYjds+AXsWAzuUug+GibOhN4TtY2xqle2HMrh9XUZpGzdwV2ha9jYZB0tnScgVDtIqQs12Oy28+gZ/m/JTtbYLW9entyDm+RLQlb9Bo5vt8rjL51u3fm06RnocJXymjGGr/ec4vW1uzCZ67kvbBWvRqQQggvixloNC7SDlKpAgzsjPFveNIsI5YWrHdzi+pzQ1R9AST50HAw3vwoDb4PwJlXvUKkg4XIblm49yrtrUul//HP+HLaaruFHcEe1IiRBO0ipqjWYhH+myMnrazOY81UmocbJi/32clPJEkKTNkJoJAy41aqE7Tws0KEqVS1FThcfpxziq7VfMPbsf3jLsZHwMCfuziPg0mcJ0Q5Sykv1PuGXlLp5Z+N+Xlm9m6iCI7wS+x1jC5bi2HsSYrppqwRVb+UVOfng6x0c+/odJpcu5Sch+3FGNMUx5F4Y/gAh2kFKVVO9Tfhut+E/aUf4y9IdxOd+y+wW6xhqvkOygd4TrHLM7mO0VYLyKRGZALwMOIDZxpiZvv6M43lFfLZiJU23zON28yXNpZD81v0wV7xE2CW36wxSqsb8kvD9fVF8k3GKV//zLf2PfcZ7EavpFH4U42iLXPVLGHY/RHfx5ccpBYCIOIB/ANcCh4AkEVlsjNlek/0tSs06r1/GtBEdaZb5H3rs/4Bpko5TwsjvdTNc8wjNdAYp5QM+T/i+vCjKXxBTL+3Cqd3f0PfQh8x1fENEmBMTezlc+mek380QGl71TpWquUuBPcaYvQAi8j4wCajRuV3W8zpOjvGT/FVMXruWVpLPyYjOnBr+NK2vvI8YLYpUPuSPO3yfXBSeF0QkxVyZt4ar1q5gUMg+SsKbEDL4bhjxU6R9fz98BaUqFAsc9Hh9CBhRkx29sCwdcZ5lfthLXONIo9SEsNydyJKIibzy5C+0KFL5hT8Svk8uiheWpVPodPGIYzEPhy6mpRSQ7u7MC46fMuM3T+tATypoich0YDpAXFxchdsczinEEEkhEbzo/BHvu0ZznBjECa9osld+ErBK26ouirJhg10I69yDead0HN+ZvkiJMEOTvQqMLMCzgqizvew8xphZwCyAxMREU9GOyuZKeMj5qwuWK+Uv/riV8PqiMMYkGmMS27ZtW371uRN/lusmfu78Gd+ZfoDoBaECKQnoJSLdRCQcmAIsrsmOZozvQ1SY47xlOnqq8jd/JHyfXBR6QahgY4wpBR4HlgE7gA+MMdtqsq/JCbE8d+sgYqOjECA2Oornbh2kg+spv/J5kY4xplREyi4KBzC3JheFDiesgpEx5gvgC1/sS0dPVXVNjKmwiLFugxA5Aey/yCZtgJN1FE5VgiWWYIkDgj+WrsaYC8sN64AX5zYEz/ELljggeGIJ9jiqdW4HRcKviogkG2MSAx0HBE8swRIHaCy1FSwxB0scEDyxNLQ4tP2XUko1EprwlVKqkagvCX9WoAPwECyxBEscoLHUVrDEHCxxQPDE0qDiqBdl+EoppWqvvtzhK6WUqiVN+Eop1UgEPOGLyAQRSReRPSLyZAXrI0Rkgb3+WxGJ91j3lL08XUTG+zmOX4nIdhHZIiKrRKSrxzqXiHxvP2rU1b6asdwnIic8PvNBj3X3ishu+3FvHcTykkccu0Qkx2Odz46LiMwVkeMisrWS9SIif7fj3CIiQz3W+fSYVCPmoDi3vYylTs7vYDm3G+15bYwJ2AOrJ24G0B0IBzYD/ctt8yjwuv18CrDAft7f3j4C6Gbvx+HHOEYDTeznj5TFYb/Or+Njch/wagXvbQXstf/G2M9j/BlLue1/htWz2h/H5RpgKLC1kvXXA0sAAS4DvvXHMalv53Ywnd/Bcm435vM60Hf458bON8aUAGVj53uaBMyzny8ExoqI2MvfN8YUG2MygT32/vwShzFmjTGmwH65EWtQOH/w5phUZjywwhiTbYw5DawAJtRhLHcC79Xi8ypljFkPZF9kk0nAfGPZCESLSEcuckxEpIuIrLHvbLeJyH/Zy1uJyAr7zmmFiMTYyyu926pAsJzbXsVSR+d3sJzbDfq8vphAJ/yKxs4vP7jIuW2MNXhVLtDay/f6Mg5P07D+1y0TKSLJIrJRRCbXMIbqxnKbnXQWikjZ6KS+PCbV2p9dBNANWO2x2JfHpSqVxXqx71AK/NoY0x/r7ukxEekPPAmsMsb0AlbZrwEmAr3sx3TgtRrEU+E2fjy3vY3Fk7/O72A5txv6eV2pejuJeaCIyFQgERjpsbirMSZLRLoDq0UkzRiT4ccwPgPeM8YUi8hDWHeJY/z4ed6YAiw0xrg8ltX1cakWY8wR4Ij9PE9EdmBdNJOAUfZm84C1wBN43G0BG0UkWkQ62vtpEILg/A62c7vendcXE9B2+CJyOfBM69atr4uPjw9YHKphS0lJKTLGRAGIyBvAWmPMeT/R7QrT9cBA4IAxJtpeLsBpY0y0iHwOzDTGfGWvWwU8YYxJLv+Zem6rupCSknLSGNO2svO6vEDf4ScBveLj40lOvuCaUconRCSkrBweuA54qtz6ZsBHwC+MMWesHG8xxhgRqdZdkfwwm9ugpk2b6rmt/EZEDtnn9gXndUUCWoZvfphQQqmaOXsKvn4ZZo+D0uLKtjqMdXORBDxrjDlXSSYiYVjJ/l1jzMf24mN2xRj23+P28mrN5gbcUtFsbkp5JWsTfPoYrHr2Ylv1o4LzujKBrrTFWBNKKOU9Y+BQMnzyMLzYD1b8ARzhkH+ssnecMsb0tB9vli20i2vmADuMMS96bL8YKGvXfC/wqcfye+zWOpcBuRcrv9dzW1VbSQGkvgOzRsO/RsPWT8BderF3bC1/Xl9MoIt0lPJeSQFsXQhJs+HIZghvBglTYfiD0L5/TfZ4JXA3kCYi39vLfgvMBD4QkWlYk5fcbq/7Aqtd9B6gALi/Nl9HqXNO7oHkufD9u1CUA236wMTnYfAUiGzps4/RhK+C38k9kDzHvhhyoV1/uOGvcMkdENG8xru1K1+lktVjK9jeAI/V+AOV8uQqhV1LrRuYvWsgJBT63QSJ0yD+KpDKTs2aqzLh2+1g5wPtAQPMMsa8LCKtgAVAPLAPuN0Yc9r+mfwy1p1QAXCfMWaTzyNXDZurFHYtsS+GtRASZl0Mwx+Erlf45WJQqk7kHYVN8yHlLTiTBS1iYfTvYOg90LyDXz/amzv8ss4pm0SkOZAiIiuwukCvMsbMtMeieBKrrbJn55QRWJ1TRvgjeNUAVXQxjPk9JNwDzdsHOjqlasYY2PeVdQOz83OrXL7HGKvYpvcEcNRNYUuVn6KdU5TfGQP7v7Yuhh2fWRdD99Fw/QvQa3ydXQxK+VxRLmx+H5LmwMl0iIyGEQ9D4gPQukedh1OtK8nunJIAfAu090jiR7GKfKDyLr+a8NX5is7AlgVWoj+x07oYLn0Ihk8LyMWglM8c2Wwl+bQPwVkAnYbCpH/AwNsgLCpgYXmd8P3YOYW4uLjqvFXVd0e3WpWwmxeA8yx0HGJdDANuhfAmgY5OqZpxFsH2RVaiP/QdhEbBoB9ZNzCdEgIdHeBlwr9Y5xRjzJGadk7BnqcxMTFR51ls6EpLYMdi627+wDcQGmkl+OEPQudhgY5OqZrLzrSaVKa+A4XZ0LonTJhpNamMiqn6/XXIm1Y6VXVOmcmFnVMeF5H3sSprL9o5RTVwOQch5U2rIvbsCYjpBtf9CYbcBU1aBTo6pWrG7YLdy60bmD2rQEKg7/VWk8ruo4K2FZk3d/jaOUVVj9sNGautYptdS61lvSdYP227j4GQgHfwVqpm8k9A6nxIfhNyD0LzjjDqSatJZYtOgY6uSt600tHOKco7BdnWz9rkuXA6E5q0gat+CcPuh+guVb9fqWBkDBzYaN3Nb/8U3E7oNhLG/xn6XA+OsEBH6DVt76ZqxxhrkKek2bD1I3AVQ9zlVkeS/jdDaESgI1SqZorz7FZkc+D4dohoaf1KTZwGbXsHOroa0YSvaqakALZ9bCX6w6kQ1hQS7rIuhg4DAx2dUjV3bJuV5LcsgJJ86HAJ3PR3q8VNeNNAR1crmvBV9ZzK+KFFQlEOtO0L1//FGtcmskWgo1OqZkqLrU5/Za3IHBFWm/nh0yB2WNBWwlaXJnxVtbJBnpLnWJWxIaHQ72brYuh6ZYO5GFQjdHq/NYxH6tuNohWZJnxVubxjdouEt+DMIWjeyR7k6V4d10bVX243ZKyy7uZ3LbNuWHpPtFuRjW7Qrcg04avzGQP7N9jj2iy2x7UZBRNnWheFjmuj6quzp6w7+eS5kLMfmraDq38Nw+5rNK3I9OpVluK8HwZ5OrHDmnTh0oesQZ7a9Ax0dEp5ZVFqFi8sS+dwTiGdoqOYcV1vJrc9bJ3X2z4GV4lVDDnuGeh7I4SGBzrkOqUJv7E7tt26my9rkdBxMNz8qlVhpePaqHpkUWoWT32cRqHTRROKuCZvFX0+XQmyH8KbW0WRw6dBu36BDjVgNOE3RufGtZkDBzbYLRLscW0aUIsE1bCVv5vPL3YSW7qfqaErudXxJS2kkB3uOJ4Pe5jf/Pp/IKJZoEMOOE34jUnuIatFQso8OHscorvCtc9Cwt0NskWCanjKknxWTiGCNQVfGKUMObOGqY6VXB6xnWITyhfuEbxTOo4U0xspEX6jyR7QhN/wud3WfJlJc6wpA42B3uNh+E+tGXcacIsE1bD8flEa7248QNnQuh04xZ2hq5jiWEs7yeGguy0znVP4wDWKbH7oE9IpOnDjzwcbTfgNVUE2bH7PSvTZGda4Nlf+lzWuTUzXQEcXFERkLnAjcNwYM9BepnM1B6FFqVm8u/EA4OaqkG3c7VjBuJAUBFjjHsLbrnGsdw/Gzfk3MFFhDmaM7xOQmIORJvyGJmuTleS3LoTSIuhymTWaX/9JOq7Nhd4CXgXmeyx7Ep2rOSh4ltHHkMc0x3rucqykW8gxTpnmzHLdyLuusRwy7c69JzoqjKYRoT+00hnfh8kJsQH8FsFFE35D4CyErWXj2myyxrUZPMWqhO0wKNDRBS1jzHp72k5POldzgHgm+JZRYeQXO+lvMvhF6ApucnxDpDhJcvfmbyW3scQ9ghLOH6UyKszBMzcP0AR/EZrw67Py49q06QMTn7eSfWROqCSCAAAf7klEQVTLQEdXX+lczQHg2aQykmKuLVnL3aEruCQkk7Mmgo9c1/COaxw7TMXFkbF6N+8VTfj1jdtldQdPmm11Dw8Jhb43WHfz8Vdrk0ofqslczaDzNV/MBR2j7CT9/NKddCw9yF2hq/iRYx0tpYB0d2f+x3kfn7iuIp8L+4QIcNdlcfxpsv6K9ZY3UxxqxVYwyD9uTROY8pbHTDu/tWfa6Rjo6BqSWs3VDDpfc2XKt7LJyinkiQ83seTDWTzvWMFVEdsoMQ6WuYfzdum1fGf6Un7uJYcIbmO0fL6GvLnDfwut2KqRyu5mvNaAZtqpR3SuZj8oa2VTluzbcZo7Hau5M3Q1HeQ0WaY1f3H+mAWu0ZwgusJ9RIU5eO7WQZrka8GbKQ61YqsGFqVmMWPhZpwu6xTPyilkxsLNAFWfsMX5HjPtbGsQM+0EIxF5D+s8biMih4Cn0bmafcqzoxQYLg/ZzlTHCsaHJBMqbta6BvN71wPccNu9zFm0g0KX69x7w0KEZpGh5BQ49Y7eR2pahq8VW1X442fbziX7Mk6X4Y+fbav8pD2+w0rym9+HkjyrhU0DmWknGBlj7qxklc7VXEuLUrN4ZvE2cgqdtOAs9zvWM9Wxkh4hRzhtmjHXNZF3XOM4YNoTGx3F7GFdkZDQ2v0iVlWqdaWtVmxdKO1QLqcLnBWuu2B5aQns/NxqbbPvS3CEw4BbrTv6zsO1ElbVO4tSs3jyoy30cO3lidAVTHJsoIkUk+ruya9LHuZz92UUY41SKXCuY9TkhFhN8H5W04SvFVsV2H0sjxdX7GLJ1qNVb5ybBZvmWZWw+cescW3GPWONa9O0jZ8jVco/cs/kkfTpP3kvZCkJoXsoMBEscl3Bu65r2Wbiz9u2rJWNJvm6U9OErxVbHg5mF/DSyl0sSs2iSXgo4we0Z9m2YxdsJ7i5ImQbvP9vSF8Cxg29rrXGtek5FkIcAYheqdpL376Zo6v/yeATn/NnyWcPnfij824+cl3DGS4sjtR284HhTbNMrdiqxLEzRbyyejcLkg4SIsKDV3fn4ZE9GPvXtedt14J8fuT4krvsMkwOtIYrHrcmF4mJD0jsStVWUXExKSsXEPX9Wwx1ptDDhLAj+hqeyhvFkqJelG9SCdrSJtC8aaWjFVvlnD5bwuvrMnhrwz5cbsOUS7vwszG9aN8i0lpvl9MPkEzudlhlmFFSwiZ3T1KGzmTYxPshLDKQX0GpGjt4IJOM5a/R59BCruQUJ6UVqd0foefERxnUNo7xqVmstXvNeoppEsbTN+nQB4GkPW2rIa/IyZyvMpn9ZSZnS0q5JSGWX4ztTVxrj16AziJuC1nP1NCVJISUlWFeyTt2Gea+m28I3BdQqoZcLjepX/6H0u9mM/Tsl3QRFzuihpF76f/S55rbaePRJ6QsoWuLm+CjCd8LRU4X87/Zx2trMzhd4GTiwA786tre9Grf/IeNsvdC8puQ+g5/Dc9mj7sTzzjv4WPX1RWWYSpVH5w6dZKtS2YRl/FvEs1BztCUrbF30OXax+jXbWCl79MWN8FJE/5FlJS6WZB8kFdW7eZ4XjEje7flv6/rw6DO9sBkbhfsXm61nd+zEiSErA5j+e/c4Xzj7k9FZZhKBTtjDNtTvyZn/esMOb2ckVLM3rBepF3yZ/peex9DI3X2qPpKE34FXG7DotQs/rZqFwezCxkeH8MrdyYwontra4P8E5D6tnVHn3vAHtfmSRh6Dzf+bSun3RW3wY+O0qEQVPA6e/Ys3y+fR/TW+Qxw7aDIhLGj9bW0Gf0o3QddHejwlA9owvdgjGHp1qP8dcUu9hzPZ2BsC/73/oGM7N3Wulc/8K09rs0icJVAt2tg/J/OG9fmdEFqpft/5uYBdfNFlKqGzN3bOLTin/Q/tpgr5QxZIR3Z1G8GfSc8RELLtoEOT/mQJnysRL9u1wn+unwXaVm59GjblNfuGsqEgR2QkrNW56ikOXAsDSJaWM0pazCujZZpqmDhdDpJXfMRoSlzGFKURBywrfmVnLpiOr0uu5FY7RPSIDX6hP9dZjZ/WZbOd/uy6RwTxV9+PJhbEmJxnNoFS16y5oUtPmOPa/MyDPqxjmuj6q1jRw6ya+lr9Nj/IZdynFNEsyl+Gj3GP8YlnboHOjzlZ4024acdyuUvy9NZt+sE7ZpH8L+TBnDH0I6E71kC8x/+YVyb/pOtyUW6XOrVuDbRUWHkFF5Yhq/l9ypQjNtN2sblFG2YxeC8dVwtpeyIGEz2sN/Sb/RdJIaFBzpEVUcaXcL3HO8mukkYT03syz0Dw4na8g688hbkH4XouBqPa/PMzQOY8eFmnO4fhgcKCxEtv1d1Ljcnm61LZ9Nh17tc4t5HHlGkdbiFTuMep1+vIYEOTwVAo0n45ce7+a8xPZked4imm5+GtV9Y49r0HAeX/t36W8MyTO10ogJtd9p3nFjzGpecWsKVUkhmaHc2DXya/uMfJLFpi0CHpwKowSf88uPdPHZ5Wx6O/o6mm38HG3ZDVCtrXJth90Orbj75TO10oupaUVEhm5e/TbO0eQxwbqWrCWVbzBhaXPMIPRJG002H2VY04ISfbY93M88e7+aXAwu5P3w1TdI+BmcBxCbC5NdhwC06ro2qt7L27WLf8n/Q5/AnjCCXw9KepF7/RZ8Jj5DQWuc6VudrcAk/r8jJ7C8zmfNVJs6SAv7QbTc/ci8lYlcKhEZZs0cNfxA6aRmmqp9cLhdb1n2CSZrN4IKNdAS2Nh3B8cum0++qyXTSJpWqEg0m4ReWWOPdvL4ug6aFWfy1/UbGFS3Dcfg0tO4JE2bC4DshquIJkpUKdtknjrBzyWvEZS4gwRwlmxakdLmXbtc9xiVxOtexqlq9T/glpW4WJB3gH6vS6VeQxJst1jHYnYTkhkDf660OUt1H6VSBql4ybjc7N60hb/0bDM5dzRXiZEf4QFKH/DcDx07l0oioQIeo6pF6m/BdbsMnqVm8uSKZq/OWsDhyDe3Cj0Foe7hmBgy7D1pqxamqn87m5ZK27E1a75hPP1cGZ00kW9reSLsxj9Gv//BAh6fqKb8kfBGZALwMOIDZxpiZvtq3MYalaUf4YulnjM5bzCLHRsLCSjFdroLEmdDvpnPj2ijla746txelZlXYdHd/eipHVv6T/sc/5zIpIDMknu/6/54BEx5keIsYn34X1fj4POGLiAP4B3AtcAhIEpHFxpjt1d2X50XRsWUktw6MIWLnR4zJ+4xXQvbjjGhKaML9MPxBpF1fX38Vpc7jq3N7UWoWT3nMCHUsJ49VH82i6+erSXBtoaNxsLXFSJpe9RC9h19Ht5AQn38X1Tj54w7/UmCPMWYvgD2h+SSgxhdFD8nirrOr+FHyelpIATkte+O++kXCBt8BETo2t6ozPjm3X1iWTqHTRVtOMzV0FXc41tBBTpNV2oaN3R+j14RHGNq+ix/CV42dPxJ+LHDQ4/UhYER1d1J2Ubwc9iqTHBsoMQ6+cI9gSeQNvPGrx7QSVgWCV+e2iEwHpgPExcVdsJPDOYUADAjZz88cn7DOfQm/cz3AGncCe++9yR9xKwUEsNLW24si1d2TdHcXFrhGcYqWiBNN9iqoGWNmAbMAEhMTTfn1naKjyMopZJ37Eq4peYlDph0AsdHa4kb5lz8KB7MAz9+jne1l5zHGzDLGJBpjEtu2vXCShU72yf+WawL/dE3iFC3PW65UAHh1bldlxvg+RIU5MIScS/ZRYQ5mjO/jmyiVqoQ/En4S0EtEuolIODAFWFzdnZRdFJ70olAB5pNze3JCLM/dOojY6CgE687+uVsH6fhLyu98XqRjjCkVkceBZVhN1+YaY7ZVdz866qQKNr46t0EH2FOBIcZcUMRY90GInAD2X2STNsDJOgqnKsESS7DEAcEfS1djTEAmZ9Vzu0aCJQ4Inlgqi6Na53ZQJPyqiEiyMSYx0HFA8MQSLHGAxlIbwRRvsMQSLHFA8MTiqzi0R4dSSjUSmvCVUqqRqC8Jf1agA/AQLLEESxygsdRGMMUbLLEESxwQPLH4JI56UYavlFKq9urLHb5SSqlaCnjCF5EJIpIuIntE5MkK1keIyAJ7/bciEu+x7il7ebqIjPdzHL8Ske0iskVEVolIV491LhH53n5UuyNODWK5T0ROeHzmgx7r7hWR3fbj3jqI5SWPOHaJSI7HOp8dFxGZKyLHRWRrJetFRP5ux7lFRIZ6rPPpMalGzHpuVz+WOjm3G+15bYwJ2AOr80oG0B0IBzYD/ctt8yjwuv18CrDAft7f3j4C6Gbvx+HHOEYDTeznj5TFYb/Or+Njch/wagXvbQXstf/G2M9j/BlLue1/htUZyR/H5RpgKLC1kvXXA0sAAS4DvvXHMdFzu/6f2435vA70Hf654WaNMSVA2XCzniYB8+znC4GxIiL28veNMcXGmExgj70/v8RhjFljjCmwX27EGkfFH7w5JpUZD6wwxmQbY04DK4AJdRjLncB7tfi8Shlj1gPZF9lkEjDfWDYC0SLSEd8fE2/puV2DWC7Cl/+Ojfa8DnTCr2i42fL9zc9tY4wpBXKB1l6+15dxeJqG9b9umUgRSRaRjSIyuYYxVDeW2+yfeAtFpGxAL18ek2rtzy4G6Aas9ljsy+NSlcpi9fUxqW08FW6j5/Z5/H1uN9rzut7OaRsoIjIVSARGeizuaozJEpHuwGoRSTPGZPgxjM+A94wxxSLyENZd4hg/fp43pgALjTEuj2V1fVxULei5XaEGdV4H+g7fm+Fmz20jIqFAS+CUl+/1ZRyIyDjgd8DNxpjisuXGmCz7715gLZBQwzi8isUYc8rj82cDw6rzPXwZi4cplPvZ6+PjUpXKYvX1MaltPBVuo+f2uc+ri3O70Z7XAW2Hb5/ku1q3bt0tPj4+YHGohi0lJcWNNfgUwCZgmDHmYuWmtVZ2bgNjsS7EJOAnxmN0TRF5DBhkjHlYRKYAtxpjbheRAcC/scqaOwGrgF7l7jJ9GUcCVh3CBGPMbo/lMUCBfbfdBvgGmGRqMD91NWLpaIw5Yj+/BXjCGHOZiLQCUrAqOKEW/47exGFv1xdYCnQzdqL09TGx9xkPfG6MGVjBuhuAx7Eqb0cAfzfGXFrj41HT2mVfPYDrhw0bZpSqlbOnKl0FlFV87gHuN3V4bmMllgzgd/ayZ7HuogEigQ/tuL4Dunu893f2+9KBiX6OYyVwDPjefiy2l18BpGG1YkkDptXBMXkO2GZ/5hqgr8d7H/DVv2NVcdivnwFmlnufT48J1q+HI4ATqxx+GvAw8LC9XoB/2HGmAYm1OR5V3uHblSbzgfaAAWYZY162/4dZAMQD+4DbjTGn7VYGL9sHtAC4zxiz6WKfkZiYaJKTky8ah1IXcLtg1zJIngP7voJfboemrS/YTERSTBCMeKhUoHlTaVsK/NoYs0lEmgMpIrICq73sKmPMTLvjwpPAE8BEoJf9GAG8Rg0mMVeqUnnHYNN8SHkLzhyC5h3hyl/oXMdKVaHKhG+s8rQj9vM8EdmB1fxnEjDK3mweVuXFE3i0GwU2iki0Z7mcUjViDOz7EpLmwM7PwV0K3UfBxJnQewI4wgIdoVJBr1rNMu3KhQTgW6C9RxI/ilXkA5W3D9WEr6qvMAc2v28V25zcBZHRMOJhGHY/tOkZ6OiUqle8Tvgi0gz4CPiFMeaMePx8NsYYEalWcx8RmQ5MB4iLi6vOW1VjcDgVkmZD2kdQWgixw2DSP2HgrRAWFejolKqXvEr4IhKGlezfNcZ8bC8+VlZUY3f1PW4v96p9qDFmFvYYz4mJiTpGs4KSAtj2sVVsc3gThDWBS26HxAeg05BAR6dUvVdlwrdb3cwBdhhjXvRYtRi4F5hp//3UY/njIvI+VmVtrpbfq4s6uRuS58L370JRLrTtCxOfh8FTILJloKNTqsHw5g7/SuBuIE1EvreX/RYr0X8gItOA/cDt9rovsJpk7sFqlnm/TyNWDYPLCelfWMU2meshJAz63QTDp0HXK7XFjVJ+4E0rna+wGv9XZGwF2xvgsVrGpRqq3CyrOeWm+ZB/FFp2gTH/A0PvgWbtAh2dUg2aDp6mfGpRahYvLEvncE4hnaKjmDG+D5MHd4S9a6xim/QvrCaWva6FxL9Br+sgxBHosJVqFDThq2qpMKEnxJ5b99THaRQ6rSFfzuYcZ9cnH5G/ch3Nzh6EJm3gip9D4v0QEx/Ab6FU46QJX3mtfELPyinkyY+2kFvopFubpvzh060UOksZKru5K3QlN4Z8S4Q4+b6gP0Num2OV0YdGBPhbKNV4acJXXlmUmsWvP9iMq9zYS0Wlbp5evI0mFDHZ8TVTw1fSP2Q/eSaKBa5RvOsayy4TR+agGwIUuVKqjCZ8VaVFqVnMWHhhsgfoLQeZ6ljJLY6vaC6F7HDH8TvnAyxyXclZrA5SsdHaUUqpYKAJX1Xpj59tw+n6IdmH42RCSBJ3ha5kRMhOSkwoJX0ns67drTy8NoRCl/vctlFhDmaM7xOIsJVS5WjCV5Uqq6A9XeAEoLOc4CeOVdzuWEsbOcM+d3ued93FwBse5frLBjISeK5V5ZW6SqnA0oSvKrQoNYsnP9pCSWkpo0M2c7djBaNCNmOAVe6hvOMax5fuQbx0x1Cu90jokxNiNcErFaQ04atzyu7os3IKaUsuDzjW8JOI1XSWkxw30bzimsz7pWM4wg+TjGhyV6r+0ITfiHm2qW8RGUp+iZNhZidPhq1gfEgS4eLia9cA/uy6ixXuYZSWO12io3QMeqXqE034jVRZyxuny9CcAiY5v2Rq6Ep6h2SRa5rwjuta3nWNJcNUfAcfFiI8c/OAOo5aKVUbmvAbqf/5dCu93ZncFbqCSY4NNJViNru7M8M5nc9cl1PEhR2kBGtS41itjFWqXtKE34gUOV0s25zJ/nXvMs+9mKEReyg04Sx2XcE7rnGkme4XvMchgtsYbXGjVAOgCb8ROJhdwJJ1XxG5ZT43udcQI/lk0JE/Ou/mI9fVnKFZhe+LCnPw3K2DNMkr1UBowm+g3G7Dl+lHSFu9gCHHFjI9ZCsuHGR3Hc/Pjozgs7yeVDTqtdjlNnpHr1TDowm/gTl9toT/bNhEyXdvMrFkOSMlm7zIdpwZ9htaXDGNts078PmT/6n0/S/dPkSTvFINlCb8BmLzgWw2rvqE+Mz3mSLJhIjhZPsrcY58hOZ9J4Djh3/qTtFRZOUUXrCP6KgwTfZKNWCa8OuxIqeLpUk7OPX1m4zO+5yHQo5yNqwluQOn03rkQ7RrdWElLMCM8X3OG+YYrPJ6bWapVMPmzSTmc4EbgePGmIH2slbAAiAe2Afcbow5bU94/jLWnLYFwH3GmE3+Cb3x2n8ynzWrl9Bq+9tMMBuIFCfHWw2h8Ko/0HTwbTQNi7zo+8vu4nXMG6UaFzEVDHl73gYi1wD5wHyPhP88kG2MmSkiTwIxxpgnROR64GdYCX8E8LIxZkRVQSQmJprk5ORafpWGzeU2rN+2j8w18xh+8hMGheyjWCI53fMW2o99DOkwKNAhBi0RSTHGJAY6DqUCzZtJzNeLSHy5xZOAUfbzecBa4Al7+Xx7IvONIhItIh2NMUd8FXBjcyq/mOXr1hG66S3Gl65htBRwslkPzlw2kxaX3kWHyBaBDlEpVU/UtAy/vUcSPwq0t5/HAgc9tjtkL7sg4YvIdGA6QFxcXA3DaJiMMaRmHidt1bv0PfgBd4bswEkYJ+LG02TMo7SJv8JuP6mUUt6rdaWtMcaIyMXLhSp+3yxgFlhFOrWNoyEoKCll5TcpFH4zmzGFyxgquZyO7MTJoU/R5qppdGrWNtAhKqXqsZom/GNlRTUi0hE4bi/PArp4bNfZXqYuYu/xM2xc/iEd97zLDWYTiHCk/dUUjXqUmL7XQUhIoENUSjUANU34i4F7gZn23089lj8uIu9jVdrmavl9xUpdbtZv3snxdbO54vRifhJygjOOGI4NeISOYx6hc7QWcymlfMubZpnvYVXQthGRQ8DTWIn+AxGZBuwHbrc3/wKrhc4erGaZ9/sh5nrtxJki1q38jBZb5zHS9Q0RUkpWzFDOXPVnWiTcQovQ8ECHqJRqoLxppXNnJavGVrCtAR6rbVANjTGGlN0HyVg5l8HHFvIjOUiBNOForzuJHfcosR36BzpEpVQjoD1t/Si/uJS161ZD8lxGFa8hUYo40rQPJ0a8QNvL76JreNNAh6iUakQ04fvBnqwTbF4+n+77FnCjpFNMOIe7TCR03ON07Dpcm1QqpQJCE76POF1uvv4uiTNfzeKq/KXcJvkcj+jMoSG/J3bUNLo1aRXoEJVSjZwm/Fo6lpPPxqXv0T79bUaZzZQSwoG2owgd/Qjt+l+rd/NKqaChCb8GjDEkb93B0TWzGHZqMZPkFKcdrcns+zhx1z5K92gdhEwpFXw04VdDXmEJG1Z+QuTmt7jC+S1h4mJvy0s5ftXztBs2mRiHHk6lVPDSDOWFXfsPsmvZLPplLWS8HCZPmpHZYypdxz9O9/a9Ax2eUkp5RRN+JUpK3Xzz1Upc3/6LywvW0ltK2N9kAAeH/5ouV99F87CoQIeolFLVogm/nMMnTrFl6Vy6ZLzHSDIoJIJ9sTfS+drH6NpNh1RXStVfmvABt9uwKfU7cta9wfDcJUyQArLCurL7kj/QY+w0+jWJDnSISilVa4064efmFZC8/B1itr9NomsLTkLZ02Y0xaMeIXbgGG1SqZRqUBplwk9P38Ghla8x6PinjJUcjoe0Y1u/X9Bz/MP0i+4Y6PCUUsovGk3CLypxkrz6I8I2vUli8bf0Ana1GEHh5dPpetlk2oU4Ah2iUkr5VYNP+FlZB0lf+jq9Di7kKo5yWlqytdv9dBv/GH079gx0eEopVWcaZMJ3u9ykblxByYZZDM1fR6w42R05iPRhT9B71E+ICYsMdIhKKVXnGlTCP306m61L/0WHXf9mmNlHPlHs6DiZTuMepVfPoYEOTymlAqpBJPz0zd9ycu0/GZy9jKulkH2h3dky8Bn6XjeNIU1aBDo8pZQKCn5J+CIyAXgZcACzjTEzff0ZRYUFfL98Ps3T5jOgdBvxJoztrcYQM/IR4geP0iaVSilVjs8Tvog4gH8A1wKHgCQRWWyM2V6T/S1KzeKFZekczimkU3QUP0sIIy5zAX2PfsplnCFLOpDc+1f0mfgwCTHtfflVlFKqQfHHHf6lwB5jzF4AEXkfmARUO+EvSs3iqY/TKHY6GROSytSzKxm5YQsGSGt2Jccu/yl9r7iJWG1SqZRSVfJHwo8FDnq8PgSMqMmOXliWTqHTxaOOxfwm7AOOmWhecU1mZdREPptxh0+CVUqpxiJglbYiMh2YDhAXF1fhNodzCgH4xHU1maYjK9zDKCUUyauzMJVSqsEI8cM+s4AuHq8728vOY4yZZYxJNMYktm3btsIddYq2hiA+QmuWuEdQav//VLZcKaWU9/yR8JOAXiLSTUTCgSnA4prsaMb4PkSFnV8+HxXmYMb4PrWPUimlGhmfF+kYY0pF5HFgGVazzLnGmG012dfkBGtuWM9WOjPG9zm3XCmllPfEGBPoGEhMTDTJycmBDkM1UCKSYozR2WtUoxcUCV9ETgD7q9isDXCyDsLxRrDEonFcqKJYuhpjKq4oUqoRCYqE7w0RSQ6Wu7RgiUXjuFAwxaJUsPFHpa1SSqkgpAlfKaUaifqU8GcFOgAPwRKLxnGhYIpFqaBSb8rwlVJK1U59usNXSilVC0GR8EVkgoiki8geEXmygvURIrLAXv+tiMR7rHvKXp4uIuP9HMevRGS7iGwRkVUi0tVjnUtEvrcfNepZXM1Y7hOREx6f+aDHuntFZLf9uNfPcbzkEcMuEcnxWOezYyIic0XkuIhsrWS9iMjf7Ti3iMhQj3U+Ox5K1WvGmIA+sHrjZgDdgXBgM9C/3DaPAq/bz6cAC+zn/e3tI4Bu9n4cfoxjNNDEfv5IWRz26/w6Pib3Aa9W8N5WwF77b4z9PMZfcZTb/mdYPav9cUyuAYYCWytZfz2wBBDgMuBbXx8Pfeijvj+C4Q7/3Pj5xpgSoGz8fE+TgHn284XAWBERe/n7xphiY0wmsMfen1/iMMasMcYU2C83Yg0M5w/eHJPKjAdWGGOyjTGngRXAhDqK407gvRp+1kUZY9YD2RfZZBIw31g2AtEi0hHfHg+l6rVgSPgVjZ9ffrCcc9sYY0qBXKC1l+/1ZRyepmHdUZaJFJFkEdkoIpNrGEN1Y7nNLr5YKCJlI5QG5JjYxVvdgNUei315TKpSWay+PB5K1WsNYhLzuiYiU4FEYKTH4q7GmCwR6Q6sFpE0Y0yGH8P4DHjPGFMsIg9h/QIa48fPq8oUYKExxuWxrK6PiVLqIoLhDt+b8fPPbSMioUBL4JSX7/VlHIjIOOB3wM3GmOKy5caYLPvvXmAtkFDDOLyKxRhzyuPzZwPDqvM9fBWHhymUK87x8TGpSmWx+vJ4KFW/BboSAetXxl6s4oCyisEB5bZ5jPMrbT+wnw/g/ErbvdS80tabOBKwKjF7lVseA0TYz9sAu7lI5aaPYuno8fwWYKP9vBWQaccUYz9v5a847O36Avuw+3X445jY+4mn8krbGzi/0vY7Xx8Pfeijvj8CXqRjKhk/X0SeBZKNMYuBOcDbIrIHq+Juiv3ebSLyAdYE6aXAY+b8IgVfx/EC0Az40Koz5oAx5magH/CGiLixfjXNNMZUe9L2asbycxG52f7e2VitdjDGZIvI/2JNRAPwrDHmYpWdtY0DrH+P940xnr34fHpMROQ9YBTQRkQOAU8DYXacrwNfYLXU2QMUAPfb63x2PJSq77SnrVJKNRLBUIavlFKqDmjCV0qpRkITvlJKNRKa8JVSqpHQhK+UUo2EJnyllGokNOErpVQjoQlfKaUaif8PUTNwdzpBAhQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 5 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.subplot(3,2,1)\n",
    "y_01 = plot_fit(CLR[:8],CL[:8])\n",
    "\n",
    "plt.subplot(3,2,2)\n",
    "y_02 =plot_fit(CLR[8:15],CL[8:15])\n",
    "\n",
    "plt.subplot(3,2,3)\n",
    "y_03 =plot_fit(CLR[15:25],CL[15:25])\n",
    "\n",
    "plt.subplot(3,2,4)\n",
    "y_04 =plot_fit(CLR[25:35],CL[25:35])\n",
    "\n",
    "plt.subplot(3,2,5)\n",
    "y_05 =plot_fit(CLR[35:45],CL[35:45])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 320,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[184.88100262 140.0086842  224.24600251 189.80924437 210.2807104 ]\n"
     ]
    }
   ],
   "source": [
    "y_hat = np.array([y_01, y_02, y_03, y_04, y_05])\n",
    "y = np.array([194, 172, 290, 186, 210])\n",
    "print(y_hat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 321,
   "metadata": {},
   "outputs": [],
   "source": [
    "def score(ground_truth, pred):\n",
    "    deltas = ground_truth - pred\n",
    "    def late(d):\n",
    "        return np.exp(-np.log(0.5)*(d/5))\n",
    "    def early(d):\n",
    "        return np.exp(np.log(0.5)*(d/20))\n",
    "    result = 0\n",
    "    for delta in deltas:\n",
    "        if delta<=0:\n",
    "            result = result + late(delta)\n",
    "        else:\n",
    "            result = result + early(delta)\n",
    "    return result/len(deltas)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 322,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "54.25959962367614\n"
     ]
    }
   ],
   "source": [
    "print(100 * score(y, y_hat))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
